using System;
using System.Collections.Generic;

using Atomic.Libraries.Mathematics;
using Atomic.Libraries.Physics.Units;
using Atomic.Thermodynamics.StateEquations;

namespace Atomic.Thermodynamics
{
	[Serializable]
	public abstract class HelmholtzPotential : IHelmholtzPotential
	{
		private IGibbsPotential gibbsPotential;

		public abstract double FreeEnergy(double temperature, double volume);

		public abstract double Pressure(double temperature, double volume);

		public abstract double PressureTemperatureDerivative(double temperature, double volume);

		public abstract double BulkModulus(double temperature, double volume);

		protected abstract IEnumerable<double> EquilibriumVolumeGuesses(double temperature);

		public double EquilibriumVolume(double temperature, double pressure)
		{
			//pressure = -0.061;

			double a = double.NegativeInfinity;
			double b = double.PositiveInfinity;

			// Search for minimum of the Gibbs free energy with a range of initial value guesses.
			foreach (double v0 in EquilibriumVolumeGuesses(temperature))
			{
				double v = v0;
				for (int i = 0; i < 100; i++)
				{
					//double f = FreeEnergy(temperature, v) + pressure * v;
					double f_v = -Pressure(temperature, v) + pressure;
					double f_vv = BulkModulus(temperature, v) / v;

					if (f_v == 0.0 && f_vv > 0.0)
					{
						// Converged exactly to solution (minimum, not maximum).
						return v;
					}

					if (-f_v > 0.0)
					{
						a = Math.Max(a, v);
					}

					if (-f_v < 0.0)
					{
						b = Math.Min(b, v);
					}

					if (!double.IsInfinity(a) && !double.IsInfinity(b))
					{
						// Found values of pressures above and below target value. Can continue using bisection.
						v = EquilibriumVolume(temperature, pressure, a, b);
						return v;
					}

					// Projected minimum by linear expansion. This solves first derivative of expansion.
					v = v - f_v / f_vv;

					if (v <= 0.0)
					{
						// Doesn't make sense to continue from here.
						break;
					}
				}

				// Try in a range around the found volume.
				v = EquilibriumVolume(temperature, pressure, 0.9 * v, 1.1 * v);
				if (!double.IsNaN(v))
				{
					return v;
				}
			}

			// Couldn't find minimum.
			return double.NaN;
		}

		public double EquilibriumVolume(double temperature, double pressure, double volumeLowerBound, double volumeUpperBound)
		{
			return RootFinder.Bisection(v => Pressure(temperature, v) - pressure, volumeLowerBound, volumeUpperBound);
		}

		/// <summary>
		/// Add several independent energy contributions (e.g. from phonons and electrons) to the same system, i.e. contained in the same volume.
		/// </summary>
		public static HelmholtzPotential Add(IEnumerable<IHelmholtzPotential> potentials)
		{
			return new AddedPotential(potentials);
		}

		/// <summary>
		/// Add several independent energy contributions (e.g. from phonons and electrons) to the same system, i.e. contained in the same volume.
		/// </summary>
		public static HelmholtzPotential Add(params IHelmholtzPotential[] potentials)
		{
			return Add((IEnumerable<IHelmholtzPotential>)potentials);
		}

		/// <summary>
		/// Extends an equation of state to the finite temperature domain, assuming temperature independent energies.
		/// </summary>
		public static HelmholtzPotential Extend(IStateEquation stateEquation)
		{
			return new StateEquationPotential(stateEquation);
		}

		/// <summary>
		/// Performs a Legendre transform to obtain the Gibbs free energy.
		/// </summary>
		public IGibbsPotential GibbsPotential
		{
			get
			{
				if (gibbsPotential == null)
				{
					gibbsPotential = new LegendreGibbsPotential(this);
				}

				return gibbsPotential;
			}
		}

		[Serializable]
		private class AddedPotential : HelmholtzPotential
		{
			private List<IHelmholtzPotential> potentials;

			public AddedPotential(IEnumerable<IHelmholtzPotential> potentials)
			{
				this.potentials = new List<IHelmholtzPotential>(potentials);
			}

			public override double FreeEnergy(double temperature, double volume)
			{
				double f = 0.0;
				foreach (IHelmholtzPotential potential in potentials)
				{
					f += potential.FreeEnergy(temperature, volume);
				}

				return f;
			}

			public override double Pressure(double temperature, double volume)
			{
				// The pressure is just minus the derivative of the free energy for each contribution. Thus we just add up.
				double p = 0.0;
				foreach (IHelmholtzPotential potential in potentials)
				{
					p += potential.Pressure(temperature, volume);
				}

				return p;
			}

			public override double PressureTemperatureDerivative(double temperature, double volume)
			{
				double p_t = 0.0;
				foreach (IHelmholtzPotential potential in potentials)
				{
					p_t += potential.PressureTemperatureDerivative(temperature, volume);
				}

				return p_t;
			}

			public override double BulkModulus(double temperature, double volume)
			{
				double b = 0.0;
				foreach (IHelmholtzPotential potential in potentials)
				{
					b += potential.BulkModulus(temperature, volume);
				}

				return b;
			}

			protected override IEnumerable<double> EquilibriumVolumeGuesses(double temperature)
			{
				foreach (IHelmholtzPotential potential in potentials)
				{
					if (potential is HelmholtzPotential)
					{
						foreach (double v in ((HelmholtzPotential)potential).EquilibriumVolumeGuesses(temperature))
						{
							yield return v;
						}
					}

					if (potential is IStateEquationPotential)
					{
						//Console.WriteLine(((IStateEquationPotential)potential).EquilibriumVolume);
						//Console.WriteLine(((IStateEquationPotential)potential).StateEquation(temperature).EquilibriumVolume);
						//throw new NotImplementedException();
						//yield return ((IStateEquationPotential)potential).EquilibriumVolume;
						yield return ((IStateEquationPotential)potential).StateEquation(temperature).EquilibriumVolume;
					}
				}
			}
		}

		[Serializable]
		private class StateEquationPotential : HelmholtzPotential
		{
			private IStateEquation stateEquation;

			public StateEquationPotential(IStateEquation stateEquation)
			{
				this.stateEquation = stateEquation;
			}

			public override double FreeEnergy(double temperature, double volume)
			{
				if (temperature < 0.0)
				{
					throw new ArgumentOutOfRangeException();
				}

				// No temperature dependence.
				return stateEquation.Energy(volume);
			}

			public override double Pressure(double temperature, double volume)
			{
				if (temperature < 0.0)
				{
					throw new ArgumentOutOfRangeException();
				}

				// No temperature dependence.
				return stateEquation.Pressure(volume);
			}

			public override double PressureTemperatureDerivative(double temperature, double volume)
			{
				if (temperature < 0.0 || volume <= 0.0)
				{
					throw new ArgumentOutOfRangeException();
				}

				// No temperature dependence.
				return 0.0;
			}

			public override double BulkModulus(double temperature, double volume)
			{
				if (temperature < 0.0)
				{
					throw new ArgumentOutOfRangeException();
				}

				// No temperature dependence.
				return stateEquation.BulkModulus(volume);
			}

			protected override IEnumerable<double> EquilibriumVolumeGuesses(double temperature)
			{
				yield return stateEquation.EquilibriumVolume;
			}
		}

		[Serializable]
		private class LegendreGibbsPotential : IGibbsPotential
		{
			private HelmholtzPotential helmholtzPotential;

			public LegendreGibbsPotential(HelmholtzPotential helmholtzPotential)
			{
				this.helmholtzPotential = helmholtzPotential;
			}

			public double FreeEnergy(double temperature, double pressure)
			{
				double volume = Volume(temperature, pressure);

				// http://en.wikipedia.org/w/index.php?title=Thermodynamic_potential&oldid=576974308#Description_and_interpretation
				return helmholtzPotential.FreeEnergy(temperature, volume) + pressure * volume;
			}

			public double Volume(double temperature, double pressure)
			{
				return helmholtzPotential.EquilibriumVolume(temperature, pressure);
			}

			public double VolumetricThermalExpansionCoefficient(double temperature, double pressure)
			{
				double volume = Volume(temperature, pressure);

				// Derived using the inverse of the Jacobian of the change-of-variables mapping (T,V)->(T,p(T,V)).
				return helmholtzPotential.PressureTemperatureDerivative(temperature, volume) / helmholtzPotential.BulkModulus(temperature, volume);
			}

			/*public double ThermalExpansionCoefficientNumerical(double temperature, double pressure)
			{
				// We have in principle all required information, as demonstated here: Here Volume depends only on the Pressure function.
				//double h = 0.001;
				//return (Volume(temperature + h, pressure) - Volume(temperature, pressure)) / h;

				double h = 0.01;

				double v = Volume(temperature, pressure);
				double vh = Volume(temperature + h, pressure);
				double alpha = 1.0 / v * (vh - v) / h;

				return alpha;
			}*/

			public IHelmholtzPotential HelmholtzPotential
			{
				get
				{
					return helmholtzPotential;
				}
			}
		}
	}
}
